function [t, y] = rk4(odefun, tspan, y0)
% rk4  使用经典四阶 Runge-Kutta 方法求解常微分方程
%
% [t, y] = rk4(odefun, tspan, y0, h)
%
% 输入参数:
%   odefun : 微分方程右侧函数句柄, 形如 f(t, y)
%   tspan  : 求解区间 [t0, tf]
%   y0     : 初始条件 (可以是标量或向量)
%   h      : 时间步长
%
% 输出参数:
%   t      : 时间向量
%   y      : 对应的数值解, 每行对应一个时间点

    t0 = tspan(1);
    tf = tspan(2);
    h = 0.001;          % 求解器的步长
    t = (t0:h:tf)';     % 构造时间向量
    N = length(t);
    
    % 初始化解数组
    y = zeros(N, length(y0));
    y(1,:) = y0;
    
    % RK4 主循环
    for i = 1:N-1
        ti = t(i);
        yi = y(i,:)';
        k1 = odefun(ti, yi);
        k2 = odefun(ti + h/2, yi + h/2 * k1);
        k3 = odefun(ti + h/2, yi + h/2 * k2);
        k4 = odefun(ti + h, yi + h * k3);
        yi_next = yi + h/6 * (k1 + 2*k2 + 2*k3 + k4);
        y(i+1,:) = yi_next';
    end
end
